library(tidyverse)
theme_set(theme_bw(base_size = 14))
library(lubridate)
rawdata <- read_tsv("data_1beep_no1st beep_annette.tsv")
rmarkdown::paged_table(rawdata)State-Space Model Analysis
1 Dataset
This dataset comes from Koval and his colleagues’ (2013) study, in which they use a novel paradigm and analytic approach to study the dynamic process of depression. In past studies, depression is associated with higher average levels of negative affect (NA) and lower average levels of positive affect (PA). However, few studies untangle the relationship between the pattern of affective fluctuations and the dynamics of depression. Therefore, in their study, they collected 99 (the final number) subjects’ depression symptoms by the Center for Epidemiologic Studies Depression Scale and (CES-D) and daily experiences of affect ratings 10 times/day for 7 days using the experience sampling method (ESM).
Affect ratings were made on:
PA = mean(happy, relaxed), and
NA = mean(sad, depressed, anxious, angry).
The scale of each affect rating is from 1 (not at all) to 100 (very much).
data <- rawdata %>%
mutate(Pos = PA,
Neg = `NA`,
Date_Time = ymd_hms(str_glue("{year}-{month}-{day} {hour}:{min}:{sec}"), tz = "CET"),
Date = as_date(Date_Time),
Time = hms::as_hms(Date_Time),
WDay = wday(Date, label = TRUE),
Subject = factor(cumsum(PpID != lag(PpID, default = 0))),
.keep = "none") %>%
group_by(Subject) %>%
mutate(Day = factor(cumsum(Date != lag(Date, default = origin)))) %>%
group_by(Subject, Date) %>%
mutate(Moment = factor(1:n())) %>%
ungroup() %>%
pivot_longer(cols = c(Pos, Neg), names_to = "Affect", values_to = "Score")
rmarkdown::paged_table(data)Nevertheless, I found some discrepancies between the data I got and the descriptions in the literature.
- The number of subjects in my data is 97;
- The sampling time points are not always 7 (days) and 10 (times). There are many missing values.
All analyses were conducted using the maximum available sample size: for analyses involving ESM data, n=95; for all other analyses, n=99.
t1 <- data %>%
group_by(Subject, Affect) %>%
summarise(num_day = length(unique(Day)),
num_obser = n(),
num_missing = sum(is.na(Score)),
num_obser_valid = num_obser - num_missing) %>%
ungroup()
rmarkdown::paged_table(t1)summary(t1) Subject Affect num_day num_obser
1 : 2 Length:194 Min. :7.000 Min. :62.00
2 : 2 Class :character 1st Qu.:7.000 1st Qu.:66.00
3 : 2 Mode :character Median :7.000 Median :66.00
4 : 2 Mean :7.031 Mean :66.75
5 : 2 3rd Qu.:7.000 3rd Qu.:67.00
6 : 2 Max. :8.000 Max. :78.00
(Other):182
num_missing num_obser_valid
Min. : 0.000 Min. :39.00
1st Qu.: 3.000 1st Qu.:58.00
Median : 6.000 Median :61.00
Mean : 6.737 Mean :60.02
3rd Qu.: 9.000 3rd Qu.:64.00
Max. :29.000 Max. :73.00
table(t1$num_day)
7 8
188 6 t2 <- data %>%
group_by(Subject, Affect, Day) %>%
summarise(num_moment = n(),
num_missing = sum(is.na(Score)),
num_valid = num_moment - num_missing) %>%
ungroup()
rmarkdown::paged_table(t2)summary(t2) Subject Affect Day num_moment
52 : 16 Length:1364 1 :194 Min. : 5.000
70 : 16 Class :character 2 :194 1st Qu.: 9.000
87 : 16 Mode :character 3 :194 Median :10.000
1 : 14 4 :194 Mean : 9.494
2 : 14 5 :194 3rd Qu.:10.000
3 : 14 6 :194 Max. :20.000
(Other):1274 (Other):200
num_missing num_valid
Min. : 0.0000 Min. : 0.000
1st Qu.: 0.0000 1st Qu.: 8.000
Median : 1.0000 Median : 9.000
Mean : 0.9582 Mean : 8.536
3rd Qu.: 1.0000 3rd Qu.:10.000
Max. :10.0000 Max. :11.000
table(t2$num_moment)
5 6 7 8 9 10 11 14 20
10 20 138 42 188 864 98 2 2 table(t2$num_valid)
0 2 3 4 5 6 7 8 9 10 11
2 4 2 30 46 89 147 208 348 448 40 2 Data Illustration
library(tsibble)subdata <- data %>%
filter(Subject %in% c(9, 20, 36, 57, 76, 85)) %>%
mutate(DateTime = as_datetime(ymd_hms(paste(as_date(as.double(Day)),
as.character(Time))))) %>%
as_tsibble(key = c(Subject, Affect),
index = DateTime)
pos_neg_color <- rev(scales::hue_pal()(2))
ggplot(subdata, aes(x = DateTime, y = Score, color = Affect)) +
geom_line() + geom_point() +
scale_y_continuous(limits = c(0, 100)) +
scale_x_datetime(breaks = as_datetime(1:7 * 86400),
labels = paste("Day", 1:7),
limits = as_datetime(c(1, 8) * 86400)) +
scale_color_manual(values = pos_neg_color) +
facet_grid(rows = "Subject") 
Some observation:
The positive affect scores generally are larger than the negative affect scores. Both scores are negatively correlated. But there are some exception cases.
The sampling time points are quite different between subjects.
The data was collected majorly in the latter half of the day.
The sampling intervals are not fixed within days and subjects.
There are a lot of missing data.
3 State-Space Model
3.1 Single Level Model
Here, I only fitted one subject specifically.
3.1.1 Model Specification
Here, I follow Schuurman & Hamaker (2019) Measurement Error Vector Autoregressive with Order 1 Model (MEVAR(1)) without a multilevel setting. It means I will fit each subject separately.
- Measurement equation
\[\begin{pmatrix}y_{1it} \\ y_{2it}\end{pmatrix} = \begin{pmatrix}\mu_{1i} \\ \mu_{2i}\end{pmatrix} + \begin{pmatrix}\theta_{1it} \\ \theta_{2it}\end{pmatrix} + \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \]
\[ \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \sim \mathcal{N} \left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \mathbf{R}_i = \begin{pmatrix} \sigma^2_{\epsilon 11i} & \sigma_{\epsilon 12i} \\ \sigma_{\epsilon 12i} & \sigma^2_{\epsilon 22i}\end{pmatrix}\right) \]
- State equation (state space model representation)
\[ \begin{pmatrix}\theta_{1it} \\ \theta_{2it}\end{pmatrix} = \begin{pmatrix} \phi_{11i} & \phi_{12i} \\ \phi_{12i} & \phi_{22i} \end{pmatrix} \begin{pmatrix}\theta_{1it-1} \\ \theta_{2it-1}\end{pmatrix} + \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \]
\[ \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \sim \mathcal{N} \left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \mathbf{Q}_i = \begin{pmatrix} \sigma^2_{\omega 11i} & \sigma_{\omega 12i} \\ \sigma_{\omega 12i} & \sigma^2_{\omega 22i}\end{pmatrix}\right) \]
where
subject index \(i = 1, \dots, N\) and occasion index: \(t = 1, \dots T\)
observation vector: \(\mathbf{y}_{it} = \begin{pmatrix} y_{1it} \\ y_{2it} \end{pmatrix}\)
latent state (within-person fluctuation): \(\boldsymbol{\theta}_{it} = \begin{pmatrix} \theta_{1it} \\ \theta_{2it}\end{pmatrix}\)
- initial latent state: \(\boldsymbol{\theta}_{i0} = \begin{pmatrix} \theta_{1i0} \\ \theta_{2i0} \end{pmatrix} \sim \mathcal{N} \left( \mathbf{m}_0 = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \mathbf{C}_0 = \begin{pmatrix} 10^3 & 0 \\ 0 & 10^3 \end{pmatrix} \right)\)
person-specific means (trait): \(\boldsymbol{\mu}_i = \begin{pmatrix}\mu_{i1} \\ \mu_{2i}\end{pmatrix}\)
measurement errors: \(\boldsymbol{\epsilon}_i = \begin{pmatrix} \epsilon_{1i} \\ \epsilon_{2i} \end{pmatrix}\)
autoregressive coefficients: \(\boldsymbol{\Phi}_i = \begin{pmatrix} \phi_{11i} & \phi_{12i} \\ \phi_{12i} & \phi_{22i} \end{pmatrix}\)
innovation: \(\boldsymbol{\omega}_i = \begin{pmatrix} \omega_{1i} \\ \omega_{2i} \end{pmatrix}\)
3.1.2 Properties
- A linear mixed-effect model expression:
We can rewrite the SSM as the LMM. If we combine the measurement equation and the state equation and extend the current time \(t\) to the origin, we have \[ \mathbf{y}_{it} = \boldsymbol{\mu}_{i} + \boldsymbol{\Phi}_i^t\boldsymbol{\theta_{i0}} + \sum_{k=0}^{t-1} \boldsymbol{\Phi}_i^k \mathbf{w}_{i,t-k} + \boldsymbol{\epsilon}_{it} \]
Let \(\mathbf{W}_{it} = \boldsymbol{\Phi}_i^t\boldsymbol{\theta_{i0}} + \sum_{k=0}^{t-1} \boldsymbol{\Phi}_i^k \mathbf{w}_{i,t-k}\), which has a mean \(0\) and a variance \(\mathbf{T}_i\) (see below). So,
\[ \mathbf{y}_{it} = \boldsymbol{\mu}_{i} + \mathbf{W}_{it} + \boldsymbol{\epsilon}_{it} \]
where \(\boldsymbol{\mu}_i\) is a fixed effect, \(\mathbf{W}_{it} \sim \mathcal{N}(\mathbf{0}, \mathbf{T}_i)\) is a Gaussian stochastic process???, and \(\epsilon_{it} \sim \mathcal{N}(0, \mathbf{R})\) is a measurement error. However, there is no random effect in this model under this case.
- The expectation of \(\mathbf{y}_{it}\)
\[\begin{align*} E[\mathbf{y}_{it}] &= E[\boldsymbol{\mu}_{i}] + E[\boldsymbol{\Phi}_i^t\boldsymbol{\theta_{i0}}] + E[\sum_{k=0}^{t-1} \boldsymbol{\Phi}_i^k \mathbf{w}_{i,t-k}] + E[\boldsymbol{\epsilon}_{it}] \\ &= \boldsymbol{\mu}_{i} + \boldsymbol{\Phi}_i^t \mathbf{m}_0 + 0 + 0 = \boldsymbol{\mu}_{i} \end{align*}\]
- The variance of \(\mathbf{y}_{it}\)
\[\begin{equation} Var(\mathbf{y}_{it}) = Var(\boldsymbol{\mu}_{i}) + Var(\boldsymbol{\theta}_{it}) + Var(\boldsymbol{\epsilon}_{it}) = 0 + Var(\boldsymbol{\theta}_{it}) + \mathbf{R}_i \end{equation}\]
\[\begin{align} Var(\boldsymbol{\theta}_{it}) = Var(\boldsymbol{\Phi}_i\mathbf{\theta}_{it-1} + \boldsymbol{w}_{it}) = \boldsymbol{\Phi}_i Var(\mathbf{\theta}_{it-1}) \boldsymbol{\Phi}_i^\top + \mathbf{Q}_i \end{align}\]
Under the stationary assumption, we let \(\mathbf{T}_i = Var(\boldsymbol{\theta}_{it})\) for all \(t\).
\[\begin{align*} &\mathbf{T}_i = \boldsymbol{\Phi}_i \mathbf{T}_i \boldsymbol{\Phi}_i^\top + \mathbf{Q}_i \\ \Rightarrow &vec(\mathbf{T}_i) = (\boldsymbol{\Phi}_i \otimes \boldsymbol{\Phi}_i) vec(\mathbf{T}_i) + vec(\mathbf{Q}_i) \\ \Rightarrow &vec(\mathbf{T}_i) = (\mathbf{I} - \boldsymbol{\Phi}_i \otimes \boldsymbol{\Phi}_i)^{-1} vec(\mathbf{Q}_i) \\ \Rightarrow &\mathbf{T}_i = mat((\mathbf{I} - \boldsymbol{\Phi}_i \otimes \boldsymbol{\Phi}_i)^{-1} vec(\mathbf{Q}_i)) \end{align*}\]
Thus,
\[ Var(\mathbf{y}_{it}) = \mathbf{T}_i + \mathbf{R}_i \]
- The covariance of \(\mathbf{y}_{it}\) and \(\mathbf{y}_{it'}\) (\(t \neq t'\)):
\[\begin{align} Cov(\mathbf{y}_{it}, \mathbf{y}_{it'}) &= Cov(\mu_i + \theta_{it} + \epsilon_{it}, \mu_i + \theta_{it'} + \epsilon_{it'}) \\ &= ????? \end{align}\]
3.1.3 Reliability (subject specific)
- Based on Schuurman & Hamaker (2019)
\[ R_k = \frac{(\mathbf{T}_i)_{kk}}{(\mathbf{T}_i)_{kk} + (\mathbf{R}_i)_{kk}} \]
where \(k = 1, 2\) for positive and negative affects, respectively.
There is a problem of representing \(Y\). In Molenberghs et al. (2007), \(\mathbf{Y}_i = (y_{i1}, \dots, y_{iT})^\top\) is a vector collecting observations for all time points. However, in our MEVAR(1) model, \(\mathbf{y}_{it} = (y_{1it}, y_{2it})^\top\) is a vector collecting the positive and negative affects at \(t\). Maybe, I should represent
\[ \mathbf{Y}_i = \begin{pmatrix} y_{1i1} & y_{2i1} \\ \vdots & \vdots \\ y_{1it} & y_{2it} \\ \vdots & \vdots \\ y_{1iT} & y_{2iT} \end{pmatrix} = \mu + ... \]
3.1.4 Bayesian Estimation by Stan
3.1.4.1 Case 1: Fit only for the positive affect
library(cmdstanr)
register_knitr_engine(override = TRUE)
library(posterior)
library(bayesplot)
color_scheme_set("red")
bayesplot_theme_set(theme_bw(base_size = 14))
# library(loo)The corresponding Stan codes for the one affect of one subject:
ssm0.stan
data {
int<lower=1> T;
vector[T] y;
real m_0;
cov_matrix[1] C_0;
}
parameters {
real mu;
real theta_0;
vector[T] theta;
real Phi;
cov_matrix[1] R;
cov_matrix[1] Q;
}
model {
theta_0 ~ normal(m_0, sqrt(C_0[1, 1]));
theta[1] ~ normal(Phi * theta_0, sqrt(Q[1, 1]));
y[1] ~ normal(mu + theta[1], sqrt(R[1, 1]));
for (t in 2:T) {
theta[t] ~ normal(Phi * theta[t-1], sqrt(Q[1, 1]));
y[t] ~ normal(mu + theta[t], sqrt(R[1, 1]));
}
}
generated quantities {
vector[T] y_hat;
y_hat = mu + theta;
}
MCMC results:
s20_pos_data <- data %>%
filter(Subject == 20, Affect == "Pos") %>%
drop_na(Score) %>%
pull(Score)
ssm0_data <- list(`T` = length(s20_pos_data),
y = s20_pos_data,
m_0 = 0,
C_0 = matrix(10^3))
ssm0 <- cmdstan_model("stan/ssm0.stan")
ssm0_fit <- ssm0$sample(data = ssm0_data,
seed = 20240207,
chains = 4,
parallel_chains = 4,
refresh = 500)
ssm0_fit$save_object(file = "stan/ssm0-fit.RDS")ssm0_fit <- readRDS("stan/ssm0-fit.RDS")
ssm0_fit$cmdstan_summary() Inference for Stan model: ssm0_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.
Warmup took (0.81, 0.85, 0.75, 0.53) seconds, 2.9 seconds total
Sampling took (0.49, 0.50, 0.69, 0.35) seconds, 2.0 seconds total
Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
lp__ -370 7.8 29 -4.0e+02 -379 -3.1e+02 13 6.6 1.3
accept_stat__ 0.83 4.9e-02 0.25 0.13 0.93 1.00 2.6e+01 1.3e+01 1.1e+00
stepsize__ 0.18 5.9e-02 0.083 0.070 0.19 0.30 2.0e+00 9.9e-01 1.5e+14
treedepth__ 4.3 3.2e-01 0.85 3.0 4.0 6.0 6.9e+00 3.4e+00 1.3e+00
n_leapfrog__ 25 6.6e+00 17 7.0 15 63 6.3e+00 3.1e+00 1.4e+00
divergent__ 0.027 nan 0.16 0.00 0.00 0.00 nan nan nan
energy__ 405 7.9e+00 29 341 414 438 1.4e+01 6.8e+00 1.3e+00
mu 55 0.13 3.3 5.0e+01 55 6.0e+01 643 318 1.0
theta_0 0.077 1.8 29 -4.8e+01 0.58 4.4e+01 244 121 1.0
theta[1] 0.79 0.22 9.7 -1.6e+01 1.1 1.6e+01 1980 978 1.0
theta[2] -4.9 0.46 9.1 -2.0e+01 -5.6 1.1e+01 398 196 1.0
theta[3] 2.7 0.24 8.5 -1.2e+01 3.1 1.6e+01 1254 619 1.0
theta[4] 4.3 0.23 8.5 -1.0e+01 4.4 1.8e+01 1407 695 1.0
theta[5] 6.8 0.36 8.9 -8.2e+00 7.3 2.2e+01 621 307 1.0
theta[6] 11 1.1 9.9 -5.7e+00 12 2.6e+01 77 38 1.1
theta[7] 4.8 0.28 8.6 -9.2e+00 5.0 1.9e+01 907 448 1.0
theta[8] 4.8 1.0 9.4 -1.1e+01 5.6 1.9e+01 85 42 1.1
theta[9] -16 1.5 11 -3.2e+01 -16 2.7e+00 55 27 1.1
theta[10] -15 0.81 10 -3.1e+01 -16 2.3e+00 163 80 1.0
theta[11] -19 1.5 11 -3.5e+01 -20 4.8e-01 60 30 1.1
theta[12] -15 0.85 10 -3.1e+01 -16 2.0e+00 143 71 1.0
theta[13] -13 1.0 10 -2.9e+01 -14 3.8e+00 96 48 1.1
theta[14] -1.1 0.20 8.3 -1.5e+01 -0.76 1.2e+01 1680 830 1.0
theta[15] 8.3 1.5 11 -1.0e+01 8.8 2.5e+01 54 27 1.1
theta[16] -11 0.70 9.3 -2.5e+01 -11 4.8e+00 177 87 1.0
theta[17] -21 2.0 13 -3.9e+01 -22 -1.1e-01 41 20 1.1
theta[18] -19 1.9 12 -3.7e+01 -20 9.3e-01 44 22 1.1
theta[19] 0.14 0.53 8.9 -1.6e+01 0.96 1.4e+01 289 143 1.0
theta[20] -3.6 0.48 8.8 -1.7e+01 -4.2 1.2e+01 344 170 1.0
theta[21] 9.6 1.1 9.6 -6.3e+00 10 2.5e+01 79 39 1.1
theta[22] 10 0.99 9.7 -6.0e+00 11 2.5e+01 96 47 1.1
theta[23] 0.24 0.30 8.7 -1.3e+01 -0.38 1.5e+01 854 422 1.0
theta[24] 9.3 1.0 9.7 -6.8e+00 9.8 2.4e+01 89 44 1.1
theta[25] 4.1 0.60 9.1 -1.1e+01 4.7 1.8e+01 230 113 1.0
theta[26] -11 1.1 10 -2.6e+01 -12 5.4e+00 82 40 1.1
theta[27] -8.2 0.27 8.6 -2.3e+01 -8.4 6.0e+00 1000 494 1.0
theta[28] -16 1.4 11 -3.3e+01 -16 2.5e+00 62 31 1.1
theta[29] -7.9 0.33 8.8 -2.2e+01 -8.2 6.7e+00 734 363 1.0
theta[30] 0.19 0.32 8.5 -1.5e+01 0.77 1.4e+01 716 353 1.0
theta[31] -3.3 0.29 8.7 -1.7e+01 -3.8 1.2e+01 910 449 1.0
theta[32] 8.2 1.2 9.7 -7.7e+00 8.7 2.3e+01 69 34 1.1
theta[33] -0.69 0.22 8.3 -1.4e+01 -1.1 1.4e+01 1432 707 1.0
theta[34] 2.9 0.19 8.6 -1.1e+01 2.7 1.7e+01 1998 987 1.0
theta[35] 21 2.9 16 -3.4e+00 20 4.4e+01 28 14 1.2
theta[36] -8.4 0.88 9.4 -2.2e+01 -9.1 7.8e+00 114 56 1.0
theta[37] -16 1.2 11 -3.2e+01 -17 2.0e+00 79 39 1.1
theta[38] -23 2.2 14 -4.3e+01 -23 -7.2e-02 38 19 1.1
theta[39] -7.9 0.26 8.7 -2.3e+01 -7.9 6.2e+00 1148 567 1.0
theta[40] 0.97 0.73 9.1 -1.5e+01 1.6 1.5e+01 156 77 1.0
theta[41] -13 1.4 11 -2.9e+01 -14 4.4e+00 56 28 1.1
theta[42] -7.4 0.40 8.9 -2.1e+01 -7.9 7.6e+00 491 243 1.0
theta[43] -1.8 0.19 8.5 -1.6e+01 -2.1 1.2e+01 2062 1018 1.0
theta[44] 10 0.86 9.5 -5.3e+00 11 2.5e+01 122 60 1.0
theta[45] 19 1.6 12 -4.4e-01 20 3.6e+01 51 25 1.1
theta[46] 14 0.79 9.8 -2.1e+00 15 2.9e+01 155 77 1.0
theta[47] 17 1.2 10 -1.8e+00 17 3.2e+01 82 41 1.1
theta[48] 14 0.86 9.9 -2.6e+00 15 3.0e+01 134 66 1.0
theta[49] 16 1.4 11 -2.3e+00 17 3.2e+01 57 28 1.1
theta[50] 8.0 0.71 9.0 -6.8e+00 8.4 2.1e+01 161 79 1.0
theta[51] -11 1.7 11 -2.8e+01 -11 7.5e+00 43 21 1.1
theta[52] -0.53 0.19 8.5 -1.5e+01 -0.37 1.3e+01 2008 992 1.00
theta[53] 9.1 1.4 10 -8.3e+00 9.4 2.5e+01 54 27 1.1
theta[54] -14 1.9 12 -3.2e+01 -13 5.2e+00 38 19 1.1
theta[55] 6.9 1.3 10 -1.0e+01 7.2 2.2e+01 62 31 1.1
theta[56] -9.0 1.5 11 -2.5e+01 -9.3 8.7e+00 48 24 1.1
theta[57] 6.3 0.62 9.2 -9.0e+00 6.7 2.1e+01 220 109 1.0
theta[58] 9.4 0.69 9.1 -5.7e+00 10 2.4e+01 173 85 1.0
theta[59] 8.0 0.32 8.8 -6.5e+00 8.4 2.3e+01 757 374 1.0
theta[60] 22 2.9 15 -5.6e-01 22 4.5e+01 26 13 1.2
theta[61] -9.8 2.0 12 -2.9e+01 -9.3 1.0e+01 36 18 1.1
theta[62] 6.1 0.59 8.7 -8.8e+00 6.7 2.0e+01 223 110 1.0
theta[63] 8.0 0.79 8.9 -6.9e+00 8.6 2.2e+01 127 63 1.0
theta[64] 0.10 0.25 8.5 -1.3e+01 -0.46 1.5e+01 1181 583 1.0
theta[65] 7.5 0.77 9.2 -7.6e+00 8.2 2.2e+01 143 71 1.0
Phi 0.36 0.016 0.24 -2.1e-02 0.36 7.6e-01 241 119 1.0
R[1,1] 178 26 120 6.7e+00 177 3.8e+02 22 11 1.2
Q[1,1] 203 24 122 2.8e+01 191 4.2e+02 26 13 1.2
y_hat[1] 56 0.16 9.3 4.0e+01 56 7.1e+01 3285 1622 1.0
y_hat[2] 50 0.55 8.8 3.6e+01 49 6.5e+01 259 128 1.0
y_hat[3] 58 0.22 8.1 4.4e+01 58 7.1e+01 1368 676 1.0
y_hat[4] 59 0.19 8.2 4.6e+01 60 7.3e+01 1792 885 1.0
y_hat[5] 62 0.32 8.6 4.7e+01 62 7.6e+01 724 358 1.0
y_hat[6] 66 1.1 9.7 5.0e+01 67 8.0e+01 78 38 1.1
y_hat[7] 60 0.24 8.2 4.6e+01 60 7.4e+01 1167 576 1.0
y_hat[8] 60 0.98 9.1 4.4e+01 61 7.3e+01 86 42 1.0
y_hat[9] 40 1.5 11 2.4e+01 39 5.7e+01 49 24 1.1
y_hat[10] 40 0.85 9.9 2.6e+01 39 5.7e+01 137 67 1.0
y_hat[11] 36 1.5 11 2.1e+01 35 5.5e+01 54 27 1.1
y_hat[12] 40 0.91 9.7 2.5e+01 39 5.7e+01 114 56 1.0
y_hat[13] 42 1.1 9.7 2.7e+01 41 5.9e+01 84 41 1.1
y_hat[14] 54 0.14 8.0 4.0e+01 54 6.7e+01 3084 1523 1.0
y_hat[15] 63 1.5 11 4.5e+01 64 7.9e+01 54 26 1.1
y_hat[16] 45 0.76 8.9 3.1e+01 44 5.9e+01 138 68 1.0
y_hat[17] 34 2.0 12 1.6e+01 33 5.5e+01 38 19 1.1
y_hat[18] 36 1.9 12 1.9e+01 35 5.6e+01 40 20 1.1
y_hat[19] 55 0.46 8.5 4.0e+01 56 6.8e+01 339 167 1.0
y_hat[20] 52 0.53 8.4 3.9e+01 51 6.6e+01 247 122 1.0
y_hat[21] 65 1.0 9.3 4.9e+01 66 7.8e+01 79 39 1.1
y_hat[22] 65 0.95 9.3 4.9e+01 66 7.9e+01 96 47 1.1
y_hat[23] 55 0.29 8.4 4.3e+01 55 7.0e+01 842 416 1.0
y_hat[24] 64 0.98 9.3 4.9e+01 65 7.8e+01 91 45 1.1
y_hat[25] 59 0.53 8.7 4.4e+01 60 7.2e+01 267 132 1.0
y_hat[26] 44 1.2 9.7 3.0e+01 43 6.0e+01 69 34 1.1
y_hat[27] 47 0.28 8.2 3.3e+01 46 6.1e+01 838 414 1.0
y_hat[28] 39 1.5 11 2.4e+01 38 5.8e+01 54 27 1.1
y_hat[29] 47 0.38 8.4 3.3e+01 47 6.1e+01 484 239 1.0
y_hat[30] 55 0.28 8.2 4.1e+01 56 6.9e+01 853 421 1.0
y_hat[31] 52 0.32 8.4 3.9e+01 51 6.6e+01 667 329 1.0
y_hat[32] 63 1.2 9.4 4.8e+01 64 7.7e+01 65 32 1.1
y_hat[33] 54 0.22 8.0 4.2e+01 54 6.9e+01 1393 688 1.0
y_hat[34] 58 0.13 8.3 4.5e+01 58 7.2e+01 4061 2006 1.0
y_hat[35] 76 2.9 15 5.2e+01 75 9.8e+01 28 14 1.2
y_hat[36] 47 0.96 9.1 3.3e+01 46 6.2e+01 90 44 1.0
y_hat[37] 39 1.2 10 2.5e+01 38 5.7e+01 70 34 1.1
y_hat[38] 32 2.2 13 1.3e+01 32 5.4e+01 35 17 1.1
y_hat[39] 47 0.23 8.3 3.4e+01 47 6.1e+01 1244 614 1.0
y_hat[40] 56 0.69 8.7 4.1e+01 57 6.9e+01 158 78 1.0
y_hat[41] 42 1.5 10 2.7e+01 41 5.9e+01 48 24 1.1
y_hat[42] 48 0.47 8.6 3.4e+01 47 6.2e+01 328 162 1.0
y_hat[43] 53 0.13 8.1 4.0e+01 53 6.7e+01 3742 1848 1.0
y_hat[44] 66 0.80 9.1 5.0e+01 66 8.0e+01 129 64 1.0
y_hat[45] 74 1.6 12 5.4e+01 75 9.0e+01 51 25 1.1
y_hat[46] 70 0.76 9.6 5.3e+01 71 8.4e+01 157 77 1.0
y_hat[47] 72 1.1 10 5.4e+01 73 8.6e+01 86 43 1.1
y_hat[48] 69 0.82 9.7 5.2e+01 71 8.4e+01 140 69 1.0
y_hat[49] 71 1.4 10 5.3e+01 72 8.6e+01 58 29 1.1
y_hat[50] 63 0.66 8.8 4.8e+01 64 7.6e+01 179 88 1.0
y_hat[51] 44 1.8 11 2.8e+01 44 6.2e+01 40 20 1.1
y_hat[52] 55 0.12 8.0 4.1e+01 55 6.8e+01 4664 2303 1.00
y_hat[53] 64 1.4 10 4.7e+01 65 7.9e+01 54 26 1.1
y_hat[54] 41 2.0 12 2.4e+01 42 6.0e+01 36 18 1.1
y_hat[55] 62 1.2 9.8 4.6e+01 62 7.6e+01 62 31 1.1
y_hat[56] 46 1.6 10 3.1e+01 46 6.3e+01 41 20 1.1
y_hat[57] 61 0.55 8.9 4.6e+01 62 7.5e+01 262 129 1.0
y_hat[58] 64 0.66 8.9 4.9e+01 65 7.9e+01 178 88 1.0
y_hat[59] 63 0.28 8.6 4.9e+01 63 7.7e+01 949 468 1.0
y_hat[60] 78 2.9 14 5.5e+01 77 9.9e+01 25 12 1.2
y_hat[61] 45 2.1 12 2.7e+01 46 6.5e+01 34 17 1.1
y_hat[62] 61 0.50 8.5 4.7e+01 62 7.4e+01 291 144 1.0
y_hat[63] 63 0.72 8.6 4.8e+01 64 7.6e+01 141 70 1.0
y_hat[64] 55 0.23 8.2 4.2e+01 54 6.9e+01 1219 602 1.0
y_hat[65] 63 0.74 9.1 4.7e+01 63 7.6e+01 150 74 1.0
Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).We can find there are sufficient ESS for each parameters and the \(\hat{R}\)’s are close to the 1.00, showing the MCMC chains could be convergent.
y_hat_pos_summary <- ssm0_fit$summary("y_hat", mean, median, quantile2)
s20_pos_predict <- subdata %>%
filter(Subject == 20, Affect == "Pos") %>%
drop_na(Score) %>%
bind_cols(y_hat_pos_summary)
ggplot(s20_pos_predict, aes(x = DateTime, y = Score)) +
geom_line(color = pos_neg_color[2]) +
geom_point(color = pos_neg_color[2]) +
geom_line(aes(y = mean), linetype = "dashed") +
geom_ribbon(aes(ymin = q5, ymax = q95), alpha = 0.25) +
scale_y_continuous(limits = c(0, 100)) +
scale_x_datetime(breaks = as_datetime(1:7 * 86400),
labels = paste("Day", 1:7),
limits = as_datetime(c(1, 8) * 86400))
The fitted line (dotted black line) with 95% CI (gray areas) indeed capture the trend of the observed data.
3.1.4.2 Case 2: Fit PA and NA simultaneously
Here, we fit the PA and NA at the same time
ssm1.stan
#include ssm-function.stan
data {
int<lower=1> T;
array[T] vector[2] y;
vector[2] m_0;
cov_matrix[2] C_0;
}
parameters {
vector[2] mu;
vector[2] theta_0;
array[T] vector[2] theta;
matrix[2, 2] Phi;
cov_matrix[2] R;
cov_matrix[2] Q;
}
model {
theta_0 ~ multi_normal(m_0, C_0);
theta[1] ~ multi_normal(Phi * theta_0, Q);
y[1] ~ multi_normal(mu + theta[1], R);
for (t in 2:T) {
theta[t] ~ multi_normal(Phi * theta[t-1], Q);
y[t] ~ multi_normal(mu + theta[t], R);
}
}
generated quantities {
array[T] vector[2] y_hat;
matrix[2, 2] Tau;
real rel_1;
real rel_2;
for (t in 1:T) {
y_hat[t] = mu + theta[t];
}
Tau = to_matrix((identity_matrix(4) - kronecker_prod(Phi, Phi)) \ to_vector(Q), 2, 2);
rel_1 = Tau[1, 1] / (Tau[1, 1] + R[1, 1]);
rel_2 = Tau[2, 2] / (Tau[2, 2] + R[2, 2]);
}
s20_data <- data %>%
filter(Subject == 20) %>%
drop_na(Score) %>%
pivot_wider(names_from = Affect, values_from = Score) %>%
select(Pos, Neg)
ssm1 <- cmdstan_model("stan/ssm1.stan")
ssm1_data <- list(`T` = 65,
y = s20_data,
m_0 = c(0, 0),
C_0 = diag(c(10^3, 10^3), 2, 2))
ssm1_fit <- ssm1$sample(data = ssm1_data,
seed = 20240207,
chains = 4,
parallel_chains = 4,
refresh = 500)
ssm1_fit$save_object(file = "stan/ssm1-fit.RDS")ssm1_fit <- readRDS("stan/ssm1-fit.RDS")
ssm1_fit$cmdstan_summary()Inference for Stan model: ssm1_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.
Warmup took (23, 25, 31, 19) seconds, 1.6 minutes total
Sampling took (18, 14, 16, 18) seconds, 1.1 minutes total
Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
lp__ -6.6e+02 2.8 30 -706 -6.7e+02 -607 111 1.7 1.0
accept_stat__ 0.76 2.3e-02 3.1e-01 0.029 0.92 1.00 176 2.7 1.0e+00
stepsize__ 0.054 9.8e-05 4.4e-03 0.047 0.054 0.060 2000 30 1.9e+13
treedepth__ 6.1 7.9e-02 9.9e-01 4.0 6.0 7.0 155 2.4 1.0e+00
n_leapfrog__ 100 3.2e+00 6.3e+01 30 63 255 401 6.1 1.0e+00
divergent__ 0.072 1.9e-02 2.6e-01 0.00 0.00 1.0 190 2.9 1.0e+00
energy__ 734 2.9e+00 3.1e+01 679 737 781 115 1.7 1.0e+00
mu[1] 5.5e+01 0.52 5.6 46 5.5e+01 64 120 1.8 1.1
mu[2] 2.5e+01 0.52 5.6 17 2.6e+01 35 115 1.8 1.1
theta_0[1] -3.5e+00 0.55 24 -44 -3.8e+00 38 1942 30 1.0
theta_0[2] -2.3e-01 0.74 21 -34 -2.0e-01 33 768 12 1.0
theta[1,1] 1.7e+00 0.62 12 -17 1.4e+00 22 371 5.7 1.0
theta[1,2] 2.4e+00 0.79 10 -15 2.6e+00 19 166 2.5 1.0
theta[2,1] -4.6e+00 0.66 11 -23 -4.3e+00 13 254 3.9 1.0
theta[2,2] -3.1e+00 0.56 9.0 -17 -3.6e+00 13 256 3.9 1.0
theta[3,1] 2.9e+00 0.60 10 -14 3.4e+00 20 283 4.3 1.0
theta[3,2] 4.0e+00 0.70 9.4 -12 4.3e+00 19 181 2.8 1.0
theta[4,1] 4.7e-01 0.58 9.6 -16 7.6e-01 15 272 4.1 1.0
theta[4,2] -2.6e+00 0.59 8.8 -16 -2.7e+00 13 223 3.4 1.0
theta[5,1] 3.8e+00 0.69 10.0 -12 3.5e+00 20 211 3.2 1.0
theta[5,2] 4.8e-01 0.77 8.9 -16 5.5e-01 15 136 2.1 1.0
theta[6,1] 5.1e+00 0.70 11 -12 4.8e+00 23 238 3.6 1.0
theta[6,2] 3.3e+00 0.79 9.5 -14 3.7e+00 19 144 2.2 1.0
theta[7,1] -5.4e-01 0.75 9.8 -16 -1.2e+00 17 171 2.6 1.0
theta[7,2] 4.1e+00 0.74 8.9 -11 4.5e+00 18 144 2.2 1.0
theta[8,1] -2.9e-01 0.73 10 -16 -5.4e-01 16 191 2.9 1.0
theta[8,2] 6.2e+00 0.84 9.2 -10.0 6.2e+00 21 120 1.8 1.0
theta[9,1] -1.4e+01 0.72 11 -33 -1.4e+01 2.2 231 3.5 1.0
theta[9,2] 1.5e+01 0.74 9.7 -1.4 1.5e+01 31 173 2.6 1.0
theta[10,1] -1.4e+01 0.85 11 -34 -1.4e+01 3.7 177 2.7 1.0
theta[10,2] 1.7e+01 0.88 9.9 1.1 1.7e+01 33 126 1.9 1.0
theta[11,1] -1.8e+01 0.71 11 -37 -1.7e+01 -0.24 257 3.9 1.0
theta[11,2] 1.8e+01 0.72 9.8 1.5 1.8e+01 34 186 2.8 1.0
theta[12,1] -1.4e+01 0.68 11 -33 -1.3e+01 2.0 257 3.9 1.0
theta[12,2] 6.8e+00 0.58 8.8 -7.0 6.6e+00 22 230 3.5 1.0
theta[13,1] -1.0e+01 0.71 11 -28 -9.4e+00 6.1 219 3.3 1.0
theta[13,2] 4.9e+00 0.60 8.9 -8.7 4.7e+00 20 217 3.3 1.0
theta[14,1] -2.5e+00 0.57 9.7 -20 -2.2e+00 12 285 4.3 1.0
theta[14,2] 4.9e+00 0.56 8.7 -9.0 4.7e+00 19 236 3.6 1.0
theta[15,1] 2.1e+00 0.77 11 -15 2.1e+00 20 204 3.1 1.0
theta[15,2] 3.7e+00 0.65 9.2 -11 3.6e+00 19 201 3.1 1.0
theta[16,1] -1.2e+01 0.64 9.8 -28 -1.1e+01 2.6 234 3.6 1.0
theta[16,2] 1.4e+01 0.67 8.6 0.20 1.4e+01 29 162 2.5 1.0
theta[17,1] -2.0e+01 0.93 13 -41 -1.9e+01 -0.62 185 2.8 1.0
theta[17,2] 2.3e+01 0.81 11 6.1 2.3e+01 41 181 2.8 1.0
theta[18,1] -2.0e+01 0.81 12 -41 -1.9e+01 -1.1 227 3.5 1.0
theta[18,2] 2.4e+01 0.77 11 6.4 2.3e+01 42 202 3.1 1.0
theta[19,1] -6.4e+00 0.87 11 -24 -6.4e+00 11 156 2.4 1.0
theta[19,2] 1.0e+01 0.58 8.8 -3.4 1.0e+01 25 226 3.4 1.0
theta[20,1] -6.3e+00 0.87 10 -23 -6.3e+00 12 135 2.1 1.0
theta[20,2] 4.0e+00 0.54 8.4 -9.5 3.8e+00 18 241 3.7 1.0
theta[21,1] 5.1e+00 0.75 11 -12 4.7e+00 24 206 3.1 1.0
theta[21,2] -3.3e+00 0.60 9.0 -17 -3.6e+00 11 225 3.4 1.0
theta[22,1] 5.6e+00 0.73 11 -11 5.0e+00 23 213 3.2 1.0
theta[22,2] -3.7e+00 0.63 9.1 -19 -3.9e+00 12 213 3.3 1.0
theta[23,1] -7.5e-01 0.66 10 -19 -7.3e-01 17 235 3.6 1.0
theta[23,2] 4.9e+00 0.64 9.1 -11 5.2e+00 19 201 3.1 1.0
theta[24,1] 4.2e+00 0.61 10 -12 3.5e+00 22 278 4.2 1.0
theta[24,2] -1.1e+00 0.59 8.8 -15 -1.5e+00 14 219 3.3 1.0
theta[25,1] 1.4e+00 0.64 10 -16 1.6e+00 19 265 4.0 1.0
theta[25,2] -5.0e+00 0.63 9.3 -20 -5.2e+00 10 219 3.3 1.0
theta[26,1] -7.5e+00 0.68 11 -26 -6.8e+00 9.1 261 4.0 1.0
theta[26,2] 9.7e+00 0.70 9.3 -5.5 9.6e+00 25 181 2.8 1.0
theta[27,1] -8.2e+00 0.52 9.6 -24 -7.7e+00 7.4 342 5.2 1.0
theta[27,2] 8.8e+00 0.50 8.3 -4.3 8.5e+00 23 269 4.1 1.0
theta[28,1] -1.2e+01 0.69 11 -31 -1.1e+01 4.4 243 3.7 1.0
theta[28,2] 7.8e+00 0.63 8.9 -6.6 7.8e+00 23 200 3.0 1.0
theta[29,1] -6.1e+00 0.67 9.8 -22 -5.9e+00 11 214 3.3 1.0
theta[29,2] 5.2e-01 0.61 8.6 -13 3.0e-01 15 200 3.1 1.0
theta[30,1] 1.6e+00 0.55 9.5 -14 2.0e+00 17 304 4.6 1.0
theta[30,2] -3.7e-01 0.55 8.3 -14 -4.1e-01 13 224 3.4 1.0
theta[31,1] -2.9e+00 0.74 11 -22 -2.5e+00 15 216 3.3 1.0
theta[31,2] -6.9e+00 0.65 9.4 -22 -7.2e+00 8.8 206 3.1 1.0
theta[32,1] 7.3e+00 0.85 12 -14 7.3e+00 28 213 3.2 1.0
theta[32,2] 6.6e+00 0.79 11 -11 6.7e+00 23 187 2.8 1.0
theta[33,1] -2.7e+00 0.93 9.8 -18 -3.1e+00 15 109 1.7 1.0
theta[33,2] -2.2e-01 0.62 8.2 -13 -3.3e-01 14 173 2.6 1.0
theta[34,1] 1.4e+00 0.60 9.1 -13 1.3e+00 17 229 3.5 1.0
theta[34,2] -2.8e+00 0.54 8.3 -16 -3.0e+00 11 230 3.5 1.0
theta[35,1] 1.3e+01 0.93 14 -7.0 1.2e+01 37 215 3.3 1.0
theta[35,2] -1.7e+00 0.78 10 -18 -1.5e+00 15 181 2.8 1.0
theta[36,1] -7.7e+00 0.68 9.9 -24 -7.1e+00 7.7 213 3.2 1.0
theta[36,2] 6.2e+00 0.61 8.6 -7.6 6.5e+00 20 199 3.0 1.0
theta[37,1] -1.3e+01 0.70 11 -33 -1.3e+01 3.1 246 3.7 1.0
theta[37,2] 8.9e+00 0.58 9.0 -5.5 8.8e+00 24 243 3.7 1.0
theta[38,1] -1.7e+01 1.2 16 -45 -1.6e+01 6.3 180 2.7 1.0
theta[38,2] 2.4e+01 1.1 14 0.93 2.4e+01 46 149 2.3 1.0
theta[39,1] -1.0e+01 0.70 10 -28 -1.0e+01 6.8 222 3.4 1.0
theta[39,2] 6.7e+00 0.61 8.7 -7.2 6.7e+00 22 201 3.1 1.0
theta[40,1] -6.7e-01 0.61 10 -17 -6.5e-01 16 277 4.2 1.0
theta[40,2] 1.5e+00 0.59 8.8 -14 1.4e+00 16 219 3.3 1.0
theta[41,1] -8.5e+00 0.76 12 -28 -7.8e+00 8.9 229 3.5 1.0
theta[41,2] 1.3e+00 0.63 9.2 -14 1.0e+00 17 217 3.3 1.0
theta[42,1] -1.7e+00 0.73 11 -20 -1.4e+00 16 224 3.4 1.0
theta[42,2] -6.0e+00 0.64 9.3 -21 -6.1e+00 9.3 212 3.2 1.0
theta[43,1] 3.0e+00 0.62 10 -15 3.2e+00 19 277 4.2 1.0
theta[43,2] -1.2e+01 0.62 9.1 -26 -1.2e+01 3.3 218 3.3 1.0
theta[44,1] 1.1e+01 0.57 9.8 -4.1 1.1e+01 27 297 4.5 1.0
theta[44,2] -1.1e+01 0.51 8.4 -25 -1.1e+01 2.6 271 4.1 1.0
theta[45,1] 1.6e+01 0.70 11 -1.4 1.5e+01 34 252 3.8 1.0
theta[45,2] -1.0e+01 0.61 9.1 -25 -1.0e+01 4.6 225 3.4 1.0
theta[46,1] 1.1e+01 0.69 10 -5.2 1.1e+01 29 225 3.4 1.0
theta[46,2] -8.2e+00 0.60 8.9 -24 -8.2e+00 6.3 221 3.4 1.0
theta[47,1] 1.3e+01 0.66 10 -3.4 1.2e+01 30 250 3.8 1.0
theta[47,2] -8.3e+00 0.53 8.7 -23 -8.2e+00 6.0 275 4.2 1.0
theta[48,1] 1.2e+01 0.81 11 -5.5 1.1e+01 32 183 2.8 1.0
theta[48,2] -1.1e+01 0.63 8.9 -25 -1.1e+01 3.7 199 3.0 1.0
theta[49,1] 1.4e+01 0.77 11 -3.6 1.3e+01 33 198 3.0 1.0
theta[49,2] -1.2e+01 0.65 9.1 -27 -1.2e+01 3.0 196 3.0 1.0
theta[50,1] 7.9e+00 0.57 9.4 -6.9 7.6e+00 24 269 4.1 1.0
theta[50,2] -7.9e+00 0.56 8.5 -22 -7.9e+00 6.5 225 3.4 1.0
theta[51,1] -4.1e+00 0.79 12 -24 -3.6e+00 14 213 3.2 1.0
theta[51,2] -2.1e+00 0.68 9.6 -18 -2.3e+00 14 200 3.0 1.0
theta[52,1] 2.5e+00 0.54 9.7 -14 2.7e+00 18 318 4.8 1.0
theta[52,2] -1.0e+01 0.59 8.8 -24 -1.0e+01 4.7 221 3.4 1.0
theta[53,1] 1.1e+01 0.60 9.9 -5.8 1.1e+01 27 274 4.2 1.0
theta[53,2] -8.9e+00 0.69 8.7 -22 -9.0e+00 5.7 159 2.4 1.0
theta[54,1] -5.7e+00 0.78 12 -28 -4.4e+00 13 251 3.8 1.0
theta[54,2] -4.7e+00 0.69 10.0 -20 -5.3e+00 13 209 3.2 1.0
theta[55,1] 9.9e+00 0.56 10 -6.7 1.0e+01 26 335 5.1 1.0
theta[55,2] -9.8e+00 0.57 8.7 -23 -1.0e+01 4.3 230 3.5 1.0
theta[56,1] -2.2e+00 0.84 13 -24 -1.3e+00 17 223 3.4 1.0
theta[56,2] -1.2e+01 0.70 10 -29 -1.2e+01 4.9 220 3.4 1.0
theta[57,1] 9.5e+00 0.48 9.9 -6.2 9.3e+00 26 413 6.3 1.0
theta[57,2] -7.7e+00 0.45 8.3 -21 -7.6e+00 5.8 346 5.3 1.0
theta[58,1] 9.1e+00 0.48 9.5 -5.9 8.3e+00 26 400 6.1 1.0
theta[58,2] -7.7e+00 0.46 8.1 -21 -7.9e+00 5.4 306 4.7 1.0
theta[59,1] 7.3e+00 0.49 9.4 -7.7 7.0e+00 24 367 5.6 1.0
theta[59,2] -1.1e+01 0.50 8.3 -25 -1.1e+01 2.4 274 4.2 1.0
theta[60,1] 1.8e+01 0.77 12 -0.31 1.7e+01 39 253 3.9 1.0
theta[60,2] -1.0e+01 0.63 9.6 -26 -1.0e+01 5.7 234 3.6 1.0
theta[61,1] -4.9e+00 0.71 12 -26 -3.8e+00 12 266 4.1 1.0
theta[61,2] -5.6e-01 0.56 9.3 -16 -7.2e-01 15 270 4.1 1.0
theta[62,1] 5.6e+00 0.51 9.3 -10 5.6e+00 21 340 5.2 1.0
theta[62,2] -6.3e+00 0.51 8.3 -20 -6.4e+00 7.4 267 4.1 1.0
theta[63,1] 5.7e+00 0.60 9.9 -10 5.3e+00 23 272 4.1 1.0
theta[63,2] -7.3e+00 0.66 8.5 -21 -7.5e+00 7.3 166 2.5 1.0
theta[64,1] 1.4e+00 0.74 12 -19 1.7e+00 21 261 4.0 1.0
theta[64,2] 9.1e+00 0.87 11 -9.3 8.8e+00 28 162 2.5 1.0
theta[65,1] 3.2e+00 0.83 12 -15 3.2e+00 22 194 3.0 1.0
theta[65,2] -7.4e+00 0.68 9.4 -23 -7.6e+00 7.8 192 2.9 1.0
Phi[1,1] -6.0e-03 0.037 0.52 -0.84 -2.0e-03 0.87 197 3.0 1.0
Phi[1,2] -6.4e-01 0.037 0.56 -1.6 -6.0e-01 0.19 229 3.5 1.0
Phi[2,1] -6.4e-02 0.038 0.54 -0.93 -7.3e-02 0.77 205 3.1 1.0
Phi[2,2] 5.8e-01 0.036 0.52 -0.26 5.7e-01 1.4 204 3.1 1.0
R[1,1] 2.5e+02 6.8 103 78 2.5e+02 420 233 3.6 1.0
R[1,2] -1.2e+02 4.8 69 -236 -1.2e+02 -14 209 3.2 1.0
R[2,1] -1.2e+02 4.8 69 -236 -1.2e+02 -14 209 3.2 1.0
R[2,2] 1.3e+02 4.3 59 33 1.2e+02 224 190 2.9 1.0
Q[1,1] 1.3e+02 7.1 95 15 1.0e+02 311 178 2.7 1.0
Q[1,2] -7.4e+01 5.2 68 -205 -6.0e+01 12 168 2.6 1.0
Q[2,1] -7.4e+01 5.2 68 -205 -6.0e+01 12 168 2.6 1.0
Q[2,2] 1.0e+02 5.2 65 18 9.3e+01 225 152 2.3 1.0
y_hat[1,1] 5.7e+01 0.31 11 40 5.7e+01 75 1188 18 1.0
y_hat[1,2] 2.8e+01 0.29 8.1 15 2.8e+01 40 762 12 1.0
y_hat[2,1] 5.1e+01 0.38 9.4 35 5.1e+01 67 594 9.1 1.0
y_hat[2,2] 2.2e+01 0.24 7.3 10 2.2e+01 35 929 14 1.0
y_hat[3,1] 5.8e+01 0.34 8.8 44 5.8e+01 73 689 11 1.0
y_hat[3,2] 2.9e+01 0.35 7.3 17 3.0e+01 41 435 6.6 1.0
y_hat[4,1] 5.6e+01 0.29 8.3 42 5.6e+01 70 808 12 1.0
y_hat[4,2] 2.3e+01 0.29 7.1 11 2.3e+01 35 595 9.1 1.0
y_hat[5,1] 5.9e+01 0.28 8.5 46 5.9e+01 74 952 15 1.0
y_hat[5,2] 2.6e+01 0.27 7.0 14 2.6e+01 37 685 10 1.0
y_hat[6,1] 6.1e+01 0.46 9.8 45 6.0e+01 77 466 7.1 1.0
y_hat[6,2] 2.9e+01 0.39 7.8 15 2.9e+01 41 408 6.2 1.0
y_hat[7,1] 5.5e+01 0.31 8.6 41 5.5e+01 69 746 11 1.0
y_hat[7,2] 3.0e+01 0.30 7.0 18 3.0e+01 41 557 8.5 1.0
y_hat[8,1] 5.5e+01 0.36 8.8 41 5.5e+01 70 594 9.1 1.0
y_hat[8,2] 3.2e+01 0.31 7.3 19 3.2e+01 43 568 8.7 1.0
y_hat[9,1] 4.1e+01 0.46 9.7 25 4.1e+01 56 441 6.7 1.0
y_hat[9,2] 4.0e+01 0.42 7.8 27 4.1e+01 53 349 5.3 1.0
y_hat[10,1] 4.1e+01 0.50 9.6 25 4.1e+01 57 377 5.7 1.0
y_hat[10,2] 4.3e+01 0.52 7.9 30 4.3e+01 55 228 3.5 1.0
y_hat[11,1] 3.8e+01 0.52 10 21 3.8e+01 54 382 5.8 1.0
y_hat[11,2] 4.3e+01 0.44 8.0 31 4.4e+01 57 334 5.1 1.0
y_hat[12,1] 4.1e+01 0.41 9.6 25 4.2e+01 57 546 8.3 1.0
y_hat[12,2] 3.2e+01 0.30 7.3 20 3.2e+01 44 588 9.0 1.0
y_hat[13,1] 4.5e+01 0.38 9.2 30 4.6e+01 60 569 8.7 1.0
y_hat[13,2] 3.0e+01 0.22 7.0 19 3.0e+01 42 1001 15 1.0
y_hat[14,1] 5.3e+01 0.30 8.4 39 5.3e+01 66 802 12 1.0
y_hat[14,2] 3.0e+01 0.21 6.7 19 3.1e+01 41 1019 16 1.0
y_hat[15,1] 5.8e+01 0.53 9.9 42 5.7e+01 75 346 5.3 1.0
y_hat[15,2] 2.9e+01 0.36 7.5 16 2.9e+01 41 441 6.7 1.0
y_hat[16,1] 4.4e+01 0.33 8.5 30 4.4e+01 58 670 10 1.0
y_hat[16,2] 4.0e+01 0.19 6.7 28 4.0e+01 51 1207 18 1.0
y_hat[17,1] 3.6e+01 0.73 11 16 3.6e+01 53 248 3.8 1.0
y_hat[17,2] 4.9e+01 0.59 9.3 34 4.8e+01 63 244 3.7 1.0
y_hat[18,1] 3.5e+01 0.65 11 16 3.6e+01 54 301 4.6 1.0
y_hat[18,2] 4.9e+01 0.58 9.5 34 4.9e+01 65 266 4.1 1.0
y_hat[19,1] 4.9e+01 0.46 9.4 34 4.9e+01 64 419 6.4 1.0
y_hat[19,2] 3.6e+01 0.24 6.9 24 3.6e+01 47 796 12 1.0
y_hat[20,1] 4.9e+01 0.38 8.4 35 4.9e+01 63 499 7.6 1.0
y_hat[20,2] 2.9e+01 0.21 6.5 19 3.0e+01 39 987 15 1.0
y_hat[21,1] 6.1e+01 0.48 9.3 46 6.0e+01 75 378 5.8 1.0
y_hat[21,2] 2.2e+01 0.35 7.3 9.9 2.2e+01 34 436 6.7 1.0
y_hat[22,1] 6.1e+01 0.54 9.7 47 6.1e+01 77 320 4.9 1.0
y_hat[22,2] 2.2e+01 0.39 7.5 9.4 2.2e+01 34 374 5.7 1.0
y_hat[23,1] 5.5e+01 0.31 8.6 40 5.5e+01 68 771 12 1.0
y_hat[23,2] 3.0e+01 0.38 7.2 18 3.0e+01 42 370 5.6 1.0
y_hat[24,1] 6.0e+01 0.46 9.3 45 5.9e+01 76 412 6.3 1.0
y_hat[24,2] 2.4e+01 0.32 7.3 13 2.4e+01 36 516 7.9 1.0
y_hat[25,1] 5.7e+01 0.50 9.6 42 5.7e+01 73 374 5.7 1.0
y_hat[25,2] 2.0e+01 0.43 8.1 7.0 2.0e+01 34 349 5.3 1.0
y_hat[26,1] 4.8e+01 0.50 9.7 32 4.8e+01 63 381 5.8 1.0
y_hat[26,2] 3.5e+01 0.43 7.5 23 3.5e+01 47 304 4.6 1.0
y_hat[27,1] 4.7e+01 0.43 8.9 32 4.8e+01 62 424 6.5 1.0
y_hat[27,2] 3.4e+01 0.24 6.8 23 3.4e+01 45 765 12 1.0
y_hat[28,1] 4.4e+01 0.50 9.5 28 4.4e+01 59 362 5.5 1.0
y_hat[28,2] 3.3e+01 0.36 7.3 21 3.3e+01 45 425 6.5 1.0
y_hat[29,1] 4.9e+01 0.38 8.7 35 5.0e+01 63 508 7.7 1.0
y_hat[29,2] 2.6e+01 0.33 7.0 15 2.5e+01 38 452 6.9 1.0
y_hat[30,1] 5.7e+01 0.23 8.0 44 5.7e+01 70 1258 19 1.0
y_hat[30,2] 2.5e+01 0.15 6.3 15 2.5e+01 35 1748 27 1.0
y_hat[31,1] 5.3e+01 0.45 9.6 36 5.3e+01 68 448 6.8 1.0
y_hat[31,2] 1.8e+01 0.42 7.9 5.6 1.9e+01 31 359 5.5 1.0
y_hat[32,1] 6.3e+01 0.65 11 45 6.2e+01 82 298 4.5 1.0
y_hat[32,2] 3.2e+01 0.53 9.0 17 3.2e+01 46 290 4.4 1.0
y_hat[33,1] 5.3e+01 0.30 8.2 39 5.3e+01 65 753 11 1.0
y_hat[33,2] 2.5e+01 0.16 6.2 15 2.5e+01 36 1488 23 1.0
y_hat[34,1] 5.7e+01 0.27 7.8 44 5.7e+01 69 834 13 1.0
y_hat[34,2] 2.3e+01 0.22 6.6 11 2.3e+01 34 928 14 1.0
y_hat[35,1] 6.9e+01 0.77 13 50 6.7e+01 92 271 4.1 1.0
y_hat[35,2] 2.4e+01 0.50 9.0 8.8 2.4e+01 38 322 4.9 1.0
y_hat[36,1] 4.8e+01 0.34 8.6 33 4.8e+01 62 646 9.8 1.0
y_hat[36,2] 3.2e+01 0.25 6.9 20 3.2e+01 43 767 12 1.0
y_hat[37,1] 4.2e+01 0.46 9.8 25 4.2e+01 58 454 6.9 1.0
y_hat[37,2] 3.4e+01 0.33 7.6 22 3.4e+01 46 536 8.2 1.0
y_hat[38,1] 3.8e+01 1.0 15 12 3.9e+01 61 214 3.3 1.0
y_hat[38,2] 4.9e+01 0.95 12 30 4.9e+01 70 168 2.6 1.0
y_hat[39,1] 4.5e+01 0.28 8.8 30 4.5e+01 59 969 15 1.0
y_hat[39,2] 3.2e+01 0.19 6.7 22 3.2e+01 44 1308 20 1.0
y_hat[40,1] 5.5e+01 0.35 8.7 40 5.5e+01 68 614 9.4 1.0
y_hat[40,2] 2.7e+01 0.29 6.9 15 2.7e+01 38 562 8.6 1.0
y_hat[41,1] 4.7e+01 0.51 10.0 30 4.8e+01 62 389 5.9 1.0
y_hat[41,2] 2.7e+01 0.35 7.4 15 2.7e+01 39 464 7.1 1.0
y_hat[42,1] 5.4e+01 0.53 9.5 38 5.4e+01 69 329 5.0 1.0
y_hat[42,2] 1.9e+01 0.39 7.7 7.7 1.9e+01 33 397 6.0 1.0
y_hat[43,1] 5.8e+01 0.51 9.4 43 5.9e+01 74 341 5.2 1.0
y_hat[43,2] 1.4e+01 0.53 8.1 1.4 1.3e+01 29 231 3.5 1.0
y_hat[44,1] 6.7e+01 0.28 8.6 53 6.6e+01 81 957 15 1.00
y_hat[44,2] 1.5e+01 0.25 7.0 3.5 1.4e+01 26 794 12 1.0
y_hat[45,1] 7.1e+01 0.48 9.9 55 7.1e+01 88 425 6.5 1.0
y_hat[45,2] 1.5e+01 0.34 7.4 3.5 1.5e+01 28 478 7.3 1.0
y_hat[46,1] 6.7e+01 0.35 8.9 52 6.7e+01 81 641 9.8 1.0
y_hat[46,2] 1.7e+01 0.29 7.1 4.9 1.7e+01 29 601 9.2 1.0
y_hat[47,1] 6.8e+01 0.38 9.1 54 6.8e+01 83 568 8.7 1.0
y_hat[47,2] 1.7e+01 0.26 7.0 5.5 1.7e+01 29 740 11 1.0
y_hat[48,1] 6.7e+01 0.52 9.5 52 6.7e+01 83 339 5.2 1.0
y_hat[48,2] 1.5e+01 0.34 7.3 2.4 1.5e+01 26 445 6.8 1.0
y_hat[49,1] 6.9e+01 0.52 9.6 54 6.9e+01 85 338 5.2 1.0
y_hat[49,2] 1.4e+01 0.35 7.3 1.6 1.3e+01 25 438 6.7 1.0
y_hat[50,1] 6.3e+01 0.29 8.2 50 6.3e+01 77 813 12 1.0
y_hat[50,2] 1.7e+01 0.22 6.8 6.6 1.7e+01 29 913 14 1.0
y_hat[51,1] 5.1e+01 0.65 11 33 5.2e+01 67 270 4.1 1.0
y_hat[51,2] 2.3e+01 0.47 8.2 11 2.3e+01 37 299 4.6 1.0
y_hat[52,1] 5.8e+01 0.41 8.8 43 5.8e+01 72 473 7.2 1.0
y_hat[52,2] 1.5e+01 0.35 7.4 3.9 1.5e+01 28 449 6.8 1.0
y_hat[53,1] 6.6e+01 0.35 8.7 52 6.6e+01 81 607 9.3 1.0
y_hat[53,2] 1.6e+01 0.31 7.0 4.8 1.6e+01 28 524 8.0 1.0
y_hat[54,1] 5.0e+01 0.69 12 30 5.1e+01 67 282 4.3 1.0
y_hat[54,2] 2.1e+01 0.49 8.7 7.2 2.0e+01 36 315 4.8 1.0
y_hat[55,1] 6.5e+01 0.27 8.8 51 6.5e+01 80 1047 16 1.0
y_hat[55,2] 1.6e+01 0.24 7.1 4.1 1.5e+01 27 890 14 1.0
y_hat[56,1] 5.3e+01 0.77 12 32 5.4e+01 72 239 3.6 1.0
y_hat[56,2] 1.3e+01 0.59 9.5 -1.6 1.3e+01 30 262 4.0 1.0
y_hat[57,1] 6.5e+01 0.48 9.1 49 6.5e+01 79 369 5.6 1.0
y_hat[57,2] 1.8e+01 0.23 6.9 6.3 1.8e+01 29 924 14 1.0
y_hat[58,1] 6.5e+01 0.38 8.9 50 6.5e+01 79 544 8.3 1.0
y_hat[58,2] 1.8e+01 0.24 6.8 6.6 1.8e+01 29 824 13 1.0
y_hat[59,1] 6.3e+01 0.39 8.8 48 6.3e+01 78 509 7.8 1.0
y_hat[59,2] 1.4e+01 0.43 7.4 2.8 1.4e+01 26 298 4.5 1.0
y_hat[60,1] 7.3e+01 0.66 11 57 7.2e+01 94 299 4.6 1.0
y_hat[60,2] 1.5e+01 0.44 8.3 1.5 1.6e+01 28 350 5.3 1.0
y_hat[61,1] 5.1e+01 0.55 11 31 5.1e+01 66 365 5.6 1.0
y_hat[61,2] 2.5e+01 0.39 7.8 12 2.5e+01 38 405 6.2 1.0
y_hat[62,1] 6.1e+01 0.37 8.7 48 6.1e+01 75 557 8.5 1.0
y_hat[62,2] 1.9e+01 0.32 7.2 7.5 1.9e+01 30 501 7.6 1.0
y_hat[63,1] 6.1e+01 0.34 8.9 47 6.1e+01 76 675 10 1.0
y_hat[63,2] 1.8e+01 0.30 7.0 6.6 1.8e+01 30 558 8.5 1.0
y_hat[64,1] 5.7e+01 0.59 11 38 5.7e+01 75 363 5.5 1.0
y_hat[64,2] 3.4e+01 0.68 9.8 19 3.4e+01 51 206 3.1 1.0
y_hat[65,1] 5.9e+01 0.55 10 42 5.9e+01 76 366 5.6 1.0
y_hat[65,2] 1.8e+01 0.39 8.2 4.0 1.8e+01 32 446 6.8 1.0
Tau[1,1] 2.6e+02 13 679 51 2.2e+02 544 2914 44 1.00
Tau[1,2] -1.9e+02 11 660 -430 -1.5e+02 -22 3649 56 1.0
Tau[2,1] -1.9e+02 11 660 -430 -1.5e+02 -22 3649 56 1.0
Tau[2,2] 2.2e+02 14 860 63 1.8e+02 463 3803 58 1.0
rel_1 5.0e-01 0.017 0.53 0.14 4.8e-01 0.86 1026 16 1.0
rel_2 6.1e-01 0.015 0.33 0.27 6.1e-01 0.92 501 7.6 1.0
Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).It seems have enough ESS for each parameters and \(\hat{R}\)’s are closed to 1.
y_hat_summary <- ssm1_fit$summary("y_hat", mean, median, quantile2) %>%
mutate(Index = str_extract(variable, "\\d+") %>% as.integer(),
Affect = str_detect(variable, ",1]") %>% ifelse("Pos", "Neg"))
s20_predict <- subdata %>%
filter(Subject == 20) %>%
drop_na(Score) %>%
mutate(Index = rep(1:(n()/2), times = 2)) %>%
left_join(y_hat_summary)
ggplot(s20_predict, aes(x = DateTime, y = Score)) +
geom_line(aes(color = Affect)) +
geom_point(aes(color = Affect)) +
scale_color_manual(values = pos_neg_color) +
geom_line(aes(y = mean, group = Affect), linetype = "dashed") +
geom_ribbon(aes(ymin = q5, ymax = q95, group = Affect), alpha = 0.25) +
scale_y_continuous(limits = c(0, 100)) +
scale_x_datetime(breaks = as_datetime(1:7 * 86400),
labels = paste("Day", 1:7),
limits = as_datetime(c(1, 8) * 86400))
The fitted lines for both PA and RA capture the observed values for Subject 20.
The posterior distributions of the reliability for positive and negative affect measurement are presented following.
ssm1_draws <- ssm1_fit$draws(format = "df") %>%
mutate(rel_diff = rel_1 - rel_2)
mcmc_areas(ssm1_draws,
pars = c("rel_1", "rel_2", "rel_diff"),
prob = 0.8,
prob_outer = 0.99)
3.2 Multilevel extension
3.2.1 Model Specification
Level 2 (between subject)
\[ \boldsymbol{\mu} = \begin{pmatrix}\mu_{1 i} \\\mu_{2 i} \end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_\mu = \begin{pmatrix} \gamma_{\mu 1} \\\gamma_{\mu 2} \end{pmatrix}, \boldsymbol{\Psi}_\mu =\begin{pmatrix} \psi_{\mu1}^2 & \psi_{\mu12} \\ \psi_{\mu12} & \psi_{\mu2}^2\end{pmatrix}\right) \]
\[ vec(\boldsymbol{\Phi}_i) = \begin{pmatrix} \phi_{11i} \\ \phi_{12i} \\ \phi_{21i} \\ \phi_{22i} \end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_{\Phi} = \begin{pmatrix}\gamma_{\phi 11} \\ \gamma_{\phi 12} \\ \gamma_{\phi 21} \\ \gamma_{\phi 22} \end{pmatrix}, \boldsymbol{\boldsymbol{\Psi}}_{\phi} = \begin{pmatrix} \psi_{\phi11}^2 & \psi_{\phi11\phi12} & \psi_{\phi11\phi21} & \psi_{\phi11\phi22} \\ \psi_{\phi11\phi12} & \psi_{\phi12}^2 & \psi_{\phi12\phi21} & \psi_{\phi12\phi22} \\ \psi_{\phi21\phi11} & \psi_{\phi21\phi12} & \psi_{\phi21}^2 & \psi_{\phi21\phi22} \\ \psi_{\phi11\phi22} & \psi_{\phi22\phi12} & \psi_{\phi22\phi21} & \psi_{\phi22}^2 \end{pmatrix} \right) \]
\[ vec(upp.tri(\mathbf{R}_i)) = \begin{pmatrix}\sigma_{\epsilon_{11 i}}^{2} \\\sigma_{\epsilon_{12 i}} \\\sigma_{\epsilon_{22 i}}^{2}\end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_{\sigma_\epsilon^2} = \begin{pmatrix} \gamma_{\sigma_{\epsilon 11}^{2}} \\\gamma_{\sigma_{\epsilon 12}} \\ \gamma_{\sigma_{\epsilon 22}^{2}} \end{pmatrix}, \boldsymbol{\Psi}_{\boldsymbol{\sigma}_{\epsilon}^{2}} = \begin{pmatrix} \psi_{\sigma_{\epsilon 11}^2} & & 0 \\ & \psi_{\sigma_{\epsilon 12}} & \\ 0 & & \psi_{\sigma_{\epsilon 22}^2}\end{pmatrix}\right) \]
\[ vec(upp.tri(\mathbf{Q}_i)) = \begin{pmatrix}\sigma_{\omega_{11 i}}^{2} \\\sigma_{\omega_{12 i}} \\\sigma_{\omega_{22 i}}^{2}\end{pmatrix} \sim \mathcal{N} \left(\boldsymbol{\gamma}_{\sigma_\omega^2} = \begin{pmatrix}\gamma_{\sigma_{\omega 11}^{2}} \\\gamma_{\sigma_{\omega 12}} \\\gamma_{\sigma_{\omega 22}^{2}} \end{pmatrix}, \boldsymbol{\Psi}_{\boldsymbol{\sigma}_{\omega}^{2}} = \begin{pmatrix} \psi_{\sigma_{\omega 11}^2} & & 0 \\ & \psi_{\sigma_{\omega 12}} & \\ 0 & & \psi_{\sigma_{\omega 22}^2}\end{pmatrix}\right) \]
However, the possible range of \(\sigma_{\epsilon 11i}^2\) is \((0, \infty)\) and the one of \(\sigma_{\epsilon 12i}\) is \((-\infty, \infty)\). Thus, the
3.2.2 Reliability
Again, according to Schuurman & Hamaker (2019), the reliability for the multilevel MEVAR(1) model have two types:
- Reliability for systematic between-subject difference
\[ R_k^B = \frac{(\Psi^2_\mu)_{kk}}{Var(\mathbf{y}_k)} = \frac{(\Psi^2_\mu)_{kk}}{(\Psi^2_\mu)_{kk} + E_{i}[(\mathbf{T}_i)_{kk}]+ (\boldsymbol{\gamma}_{\sigma_\epsilon^2})_k} \]
- Reliability for within-subject fluctuations
\[ R_{ik}^W = \frac{(\mathbf{T}_i)_{kk}}{Var(y_{ik})} = \frac{(\mathbf{T}_i)_{kk}}{(\mathbf{T}_i)_{kk} + (\mathbf{R}_i)_{kk}} \]
3.2.3 Bayesian Estimation
multilevel-ssm.stan
#include ssm-function.stan
data {
int<lower=1> N; // number of subjects
array[N] int<lower=1> T; // number of observation for each subject
int<lower=1> max_T; // maximum number of observation
int<lower=1> P; // number of affects
array[N] matrix[P, max_T] y; // observations
array[N] vector[P] m_0; // prior mean of the intial state
array[N] cov_matrix[P] C_0; // prior covariance of the intial state
}
parameters {
array[N] vector[P] mu; // ground mean/trane
array[N] vector[P] theta_0; // initial latent state
array[N] matrix[P, max_T] theta; // latent states
array[N] matrix[P, P] Phi; // autoregressive parameters
array[N] cov_matrix[P] R; // covariance of the measurment error
array[N] cov_matrix[P] Q; // covariance of the innovation noise
vector[P] gamma_mu; // prior mean of the ground mean
cov_matrix[P] Psi_mu; // prior covariance of the ground mean
vector[P * P] gamma_Phi; // prior mean of the autoregressive parameters
cov_matrix[P * P] Psi_Phi; // prior covariance of the autoregressive parameters
vector[P * (P + 1) / 2] gamma_R; // prior mean of the covariance of the measurement error
vector[P * (P + 1) / 2] diag_Psi_R; // prior covariance of the covariance of the measurement error
vector[P * (P + 1) / 2] gamma_Q; // prior mean of the covariance of the innovation noise
vector[P * (P + 1) / 2] diag_Psi_Q; // prior covariance of the covariance of innovation noise
}
model {
// level 1 (within subject)
for (n in 1:N) {
// when t = 0
theta_0[n] ~ multi_normal(m_0[n], C_0[n]);
// when t = 1
theta[n][, 1] ~ multi_normal(Phi[n] * theta_0[n], Q[n]);
y[n][, 1] ~ multi_normal(mu[n] + theta[n][, 1], R[n]);
// when t = 2, ..., T
for (t in 2:T[n]) {
theta[n][, t] ~ multi_normal(Phi[n] * theta[n][, t - 1], Q[n]);
y[n][, t] ~ multi_normal(mu[n] + theta[n][, t], R[n]);
}
}
// level 2 (between subject)
for (n in 1:N) {
mu[n] ~ multi_normal(gamma_mu, Psi_mu);
to_vector(Phi[n]) ~ multi_normal(gamma_Phi, Psi_Phi);
to_vector_lower_tri(R[n]) ~ normal(gamma_R, sqrt(diag_Psi_R));
to_vector_lower_tri(Q[n]) ~ normal(gamma_Q, sqrt(diag_Psi_Q));
}
// the (hyper)priors of parameters are set as the Stan default values
}
generated quantities {
array[N] matrix[P, max_T] y_hat;
array[N] matrix[P, P] Tau;
array[N] vector[P] rel_W;
vector[P] rel_B;
for (n in 1:N) {
// prediction
for (t in 1:T[n]) {
y_hat[n][, t] = mu[n] + theta[n][, t];
}
// within-subject reliability
Tau[n] = to_matrix((identity_matrix(P * P) - kronecker_prod(Phi[n], Phi[n])) \ to_vector(Q[n]), P, P);
for (p in 1:P) {
rel_W[n, p] = Tau[n, p, p] / (Tau[n, p, p] + R[n, p, p]);
}
}
// between-subject reliability
for (p in 1:P) {
rel_B[p] = Psi_mu[p, p] / (Psi_mu[p, p] + mean(Tau[, p, p]) + gamma_R[index_of_diag_lower_tri(p, P)]);
}
}
data_list <- data %>%
group_by(Subject) %>%
pivot_wider(names_from = Affect, values_from = Score) %>%
select(Pos, Neg) %>%
drop_na(Pos, Neg) %>%
nest() %>% ungroup() %>%
mutate(`T` = map_dbl(data, nrow),
max_T = max(`T`),
data_padding = pmap(list(data, `T`, max_T),
\(x, y, z) {
bind_rows(x,
tibble(Pos = rep(Inf, z - y),
Neg = rep(Inf, z - y))) %>%
t()
}))
mssm <- cmdstan_model("stan/multilevel-ssm.stan")
mssm_data <- lst(N = nrow(data_list),
`T` = map(data_list$data, nrow),
max_T = max(data_list$`T`),
P = 2,
y = data_list$data_padding,
m_0 = rep(list(c(0, 0)), N),
C_0 = rep(list(diag(c(10^3, 10^3), 2, 2)), N))
mssm_fit <- mssm$sample(data = mssm_data,
seed = 20240221,
chains = 4,
parallel_chains = 4,
iter_sampling = 3000,
refresh = 1000,
show_messages = FALSE)
mssm_fit$save_object(file = "stan/multilevel-ssm-fit.RDS")It took me more than 15 hours to finish the sampling …
mssm_fit <- readRDS("stan/multilevel-ssm-fit.RDS")
mssm_fit_summary <- readRDS("stan/multilevel-ssm-summary.RDS")
skimr::skim(data.frame(Rhat = as.double(mssm_fit_summary$rhat),
ESS = as.double(mssm_fit_summary$ess_bulk)))| Name | data.frame(Rhat = as.doub… |
| Number of rows | 30499 |
| Number of columns | 2 |
| _______________________ | |
| Column type frequency: | |
| numeric | 2 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Rhat | 2522 | 0.92 | 1.27 | 0.40 | 1.01 | 1.07 | 1.12 | 1.24 | 3.93 | ▇▁▁▁▁ |
| ESS | 2522 | 0.92 | 39.90 | 45.39 | 4.30 | 12.25 | 23.14 | 48.24 | 562.51 | ▇▁▁▁▁ |
However, the most of the ESS are super low. And \(\hat{R}\)’s are deviated form 1, meaning the MCMC samples are not converged yet. Thus, the posterior results are not reliable.
y_hat_summary <- mssm_fit$summary("y_hat", mean, median, quantile2) %>%
mutate(Indices = str_extract_all(variable, "\\d+"),
Subject = map_dbl(Indices, \(x) as.integer(x[1])) %>%
factor(levels = levels(data$Subject)),
Affect = map_chr(Indices, \(x) ifelse(x[2] == 1, "Pos", "Neg")),
Time_Index = map_dbl(Indices, \(x) as.integer(x[3])))
data_predict <- data %>%
mutate(DateTime = as_datetime(ymd_hms(paste(as_date(as.double(Day)),
as.character(Time))))) %>%
as_tsibble(key = c(Subject, Affect),
index = DateTime) %>%
pivot_wider(names_from = Affect, values_from = Score) %>%
select(Pos, Neg) %>%
drop_na(Pos, Neg) %>%
group_by(Subject) %>%
mutate(Time_Index = 1:n()) %>%
ungroup() %>%
pivot_longer(c("Pos", "Neg"), names_to = "Affect", values_to = "Score") %>%
left_join(y_hat_summary)
for (i in 1:10) {
g <- data_predict %>%
filter(Subject %in% (10 * (i - 1) + 1):(10 * i)) %>%
ggplot(aes(x = DateTime, y = Score)) +
geom_line(aes(color = Affect)) +
geom_point(aes(color = Affect)) +
scale_color_manual(values = pos_neg_color) +
geom_line(aes(y = mean, group = Affect), linetype = "dashed") +
geom_ribbon(aes(ymin = q5, ymax = q95, group = Affect), alpha = 0.25) +
geom_hline(yintercept = c(0, 100), color = "forestgreen") +
scale_y_continuous(limits = c(-20, 120)) +
scale_x_datetime(breaks = as_datetime(1:7 * 86400),
labels = paste("Day", 1:7),
limits = as_datetime(c(1, 8) * 86400)) +
facet_grid(Subject ~ .)
ggsave(filename = str_glue("plots/mssm2/predict_{from}-{to}.png",
from = 10 * (i - 1) + 1, to = 10 * i),
plot = g, width = 7, height = 14)
}The fitting recovery results for Subjects 1-20. See the other subject results from this link.
3.2.4 Within- and Between-Subject Reliability
rel_W_draws <- mssm_fit$draws(variables = "rel_W", format = "df") %>%
select(-.chain, -.iteration, -.draw) %>%
pivot_longer(cols = everything(),
names_to = "variable", values_to = "value") %>%
mutate(Indices = str_extract_all(variable, "\\d+"),
Subject = map_dbl(Indices, \(x) as.integer(x[1])) %>%
factor(levels = levels(data$Subject)),
Affect = map_chr(Indices, \(x) ifelse(x[2] == 1, "Pos", "Neg")))
g_rel_W <- ggplot(rel_W_draws, aes(x = 1, y = value, fill = Affect)) +
geom_split_violin() +
scale_x_continuous(name = NULL, labels = NULL, breaks = NULL) +
scale_y_continuous(limits = c(0, 1)) +
scale_fill_manual(values = pos_neg_color) +
facet_wrap(~ Subject, ncol = 10)
ggsave("plots/rel-W.png", g_rel_W, width = 14, height = 14)
rel_B_draws <- mssm_fit$draws(variables = "rel_B", format = "df") %>%
mutate(rel_B_diff = `rel_B[1]` - `rel_B[2]`)
g_rel_B <- mcmc_areas(rel_B_draws,
prob = 0.8,
prob_outer = 0.99) +
coord_cartesian(xlim = c(-1.5, 1.5))
ggsave("plots/rel-B.png", g_rel_B)
3.3 Re-parameterize the model
Refer from: https://mc-stan.org/docs/stan-users-guide/multivariate-hierarchical-priors.html
\[ \mathbf{R} = diag.matrix(\boldsymbol{\tau}_R) \times \boldsymbol{\Omega} \times diag.matrix(\boldsymbol{\tau}_R) \]
where \(\boldsymbol{\tau}_R\) is the vector of coefficient scale, \(\boldsymbol{\Omega}\) is the correlation matrix, which can be decomposed by a product of the lower triangle matrix \(\boldsymbol{\Omega_L}\)through Cholesky factorization, \(\boldsymbol{\Omega}_R = \boldsymbol{\Omega_{RL} \Omega_{RL}'}\), for optimization.
\[\begin{align*} \boldsymbol{\tau}_R &\sim Cauchy(loc_{\tau_R}, scale_{\tau_R}) \\ \boldsymbol{\Omega}_{RL} &\sim LKJ.corr.cholesky(\eta_{\Omega_R}) \end{align*}\]
The same setting is applied for \(\mathbf{Q}\).
multilevel-ssm.stan
#include ssm-function.stan
data {
int<lower=1> N; // number of subjects
array[N] int<lower=1> T; // number of observation for each subject
int<lower=1> max_T; // maximum number of observation
int<lower=1> P; // number of affects
array[N] matrix[P, max_T] y; // observations
array[N] vector[P] m_0; // prior mean of the intial state
array[N] cov_matrix[P] C_0; // prior covariance of the intial state
}
parameters {
array[N] vector[P] mu; // ground mean/trane
array[N] vector[P] theta_0; // initial latent state
array[N] matrix[P, max_T] theta; // latent states
array[N] matrix[P, P] Phi; // autoregressive parameters
array[N] cholesky_factor_corr[P] L_Omega_R;
array[N] cholesky_factor_corr[P] L_Omega_Q;
array[N] vector<lower=0>[P] tau_R;
array[N] vector<lower=0>[P] tau_Q;
vector[P] gamma_mu; // prior mean of the ground mean
cov_matrix[P] Psi_mu; // prior covariance of the ground mean
vector[P * P] gamma_Phi; // prior mean of the autoregressive parameters
cov_matrix[P * P] Psi_Phi; // prior covariance of the autoregressive parameters
vector<lower=0>[P] loc_tau_R; // location of the measurement error
vector<lower=0>[P] loc_tau_Q; // location of the innovation noise
vector<lower=0>[P] scale_tau_R; // scale of the measurement error
vector<lower=0>[P] scale_tau_Q; // scale of the innovation noise
real<lower=0> eta_Omega_R; // shape of the LKJ prior for the measurement error
real<lower=0> eta_Omega_Q; // shape of the LKJ prior for the innovation noise
}
transformed parameters{
array[N] cov_matrix[P] R; // covariance of the measurment error
array[N] cov_matrix[P] Q; // covariance of the innovation noise
for (n in 1:N) {
R[n] = diag_pre_multiply(tau_R[n], L_Omega_R[n]) * diag_pre_multiply(tau_R[n], L_Omega_R[n])';
Q[n] = diag_pre_multiply(tau_Q[n], L_Omega_Q[n]) * diag_pre_multiply(tau_Q[n], L_Omega_Q[n])';
}
}
model {
// level 1 (within subject)
for (n in 1:N) {
// when t = 0
theta_0[n] ~ multi_normal(m_0[n], C_0[n]);
// when t = 1
theta[n][, 1] ~ multi_normal(Phi[n] * theta_0[n], Q[n]);
y[n][, 1] ~ multi_normal(mu[n] + theta[n][, 1], R[n]);
// when t = 2, ..., T
for (t in 2:T[n]) {
theta[n][, t] ~ multi_normal(Phi[n] * theta[n][, t - 1], Q[n]);
y[n][, t] ~ multi_normal(mu[n] + theta[n][, t], R[n]);
}
}
// level 2 (between subject)
for (n in 1:N) {
mu[n] ~ multi_normal(gamma_mu, Psi_mu);
to_vector(Phi[n]) ~ multi_normal(gamma_Phi, Psi_Phi);
tau_R[n] ~ cauchy(loc_tau_R, scale_tau_R);
tau_Q[n] ~ cauchy(loc_tau_Q, scale_tau_Q);
L_Omega_R[n] ~ lkj_corr_cholesky(eta_Omega_R);
L_Omega_Q[n] ~ lkj_corr_cholesky(eta_Omega_Q);
}
// the (hyper)priors of parameters are set as the Stan default values
}
generated quantities {
array[N] matrix[P, max_T] y_hat;
array[N] matrix[P, P] Tau;
array[N] vector[P] rel_W;
vector[P] rel_B;
for (n in 1:N) {
// prediction
for (t in 1:T[n]) {
y_hat[n][, t] = mu[n] + theta[n][, t];
}
// within-subject reliability
Tau[n] = to_matrix((identity_matrix(P * P) - kronecker_prod(Phi[n], Phi[n])) \ to_vector(Q[n]), P, P);
for (p in 1:P) {
rel_W[n, p] = Tau[n, p, p] / (Tau[n, p, p] + R[n, p, p]);
}
}
// between-subject reliability
for (p in 1:P) {
rel_B[p] = Psi_mu[p, p] / (Psi_mu[p, p] + mean(Tau[, p, p]) + loc_tau_R[p]^2);
}
}
mssm2 <- cmdstan_model("stan/multilevel-ssm2.stan")
data_list <- data %>%
#filter(Subject %in% c(1:9, 87)) %>%
group_by(Subject) %>%
pivot_wider(names_from = Affect, values_from = Score) %>%
select(Pos, Neg) %>%
drop_na(Pos, Neg) %>%
nest() %>% ungroup() %>%
mutate(`T` = map_dbl(data, nrow),
max_T = max(`T`),
data_padding = pmap(list(data, `T`, max_T),
\(x, y, z) {
bind_rows(x,
tibble(Pos = rep(Inf, z - y),
Neg = rep(Inf, z - y))) %>%
t()
}))
mssm_data <- lst(N = nrow(data_list),
`T` = map(data_list$data, nrow),
max_T = max(data_list$`T`),
P = 2,
y = data_list$data_padding,
m_0 = rep(list(c(0, 0)), N),
C_0 = rep(list(diag(c(10^3, 10^3), 2, 2)), N))
mssm2_fit <- mssm2$sample(data = mssm_data,
seed = 20240221,
chains = 8,
parallel_chains = 8,
iter_sampling = 1500,
refresh = 1000,
show_messages = FALSE)
mssm2_fit$save_object(file = "stan/multilevel-ssm2-fit.RDS")
saveRDS(mssm2_fit$summary(), "stan/multilevel-ssm2-summary.RDS")| Name | data.frame(Rhat = as.doub… |
| Number of rows | 31661 |
| Number of columns | 2 |
| _______________________ | |
| Column type frequency: | |
| numeric | 2 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Rhat | 3002 | 0.91 | 1.55 | 0.46 | 1.00 | 1.31 | 1.4 | 1.56 | 5.12 | ▇▁▁▁▁ |
| ESS | 3002 | 0.91 | 31.39 | 172.43 | 8.37 | 14.76 | 19.0 | 25.70 | 12032.13 | ▇▁▁▁▁ |
rel_W_draws_2 <- mssm2_fit$draws(variables = "rel_W", format = "df") %>%
select(-.chain, -.iteration, -.draw) %>%
pivot_longer(cols = everything(),
names_to = "variable", values_to = "value") %>%
mutate(Indices = str_extract_all(variable, "\\d+"),
Subject = map_dbl(Indices, \(x) as.integer(x[1])) %>%
factor(levels = levels(data$Subject)),
Affect = map_chr(Indices, \(x) ifelse(x[2] == 1, "Pos", "Neg")))
g_rel_W_2 <- ggplot(rel_W_draws_2, aes(x = 1, y = value, fill = Affect)) +
geom_split_violin() +
scale_x_continuous(name = NULL, labels = NULL, breaks = NULL) +
scale_y_continuous(limits = c(0, 1)) +
scale_fill_manual(values = pos_neg_color) +
facet_wrap(~ Subject, ncol = 10)
ggsave("plots/rel-W-2.png", g_rel_W_2, width = 14, height = 14)
rel_B_draws_2 <- mssm2_fit$draws(variables = "rel_B", format = "df") %>%
mutate(rel_B_diff = `rel_B[1]` - `rel_B[2]`)
g_rel_B_2 <- mcmc_areas(rel_B_draws_2,
prob = 0.8,
prob_outer = 0.99) +
coord_cartesian(xlim = c(-1.5, 1.5))
ggsave("plots/rel-B-2.png", g_rel_B_2)
However, the results are worse than the original one. The main reason is the MCMC samples had not converged yet.
References
Koval, P., Pe, M. L., Meers, K., & Kuppens, P. (2013). Affect dynamics in relation to depressive symptoms: Variable, unstable or inert? Emotion, 13(6), 1132–1141. https://doi.org/10.1037/a0033579
Molenberghs, G., Laenen, A., & Vangeneugden, T. (2007). Estimating reliability and generalizability from hierarchical biomedical data. Journal of Biopharmaceutical Statistics, 17(4), 595–627. https://doi.org/10.1080/10543400701329448
Schuurman, N. K., & Hamaker, E. L. (2019). Measurement error and person-specific reliability in multilevel autoregressive modeling. Psychological Methods, 24(1), 70–91. https://doi.org/10.1037/met0000188
Vangeneugden, T., Laenen, A., Geys, H., Renard, D., & Molenberghs, G. (2005). Applying concepts of generalizability theory on clinical trial data to investigate sources of variation and their impact on reliability. Biometrics, 61(1), 295–304. https://doi.org/10.1111/j.0006-341X.2005.031040.x